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Abstract 

This study examines the simulation of quantum algorithms on a classical 
computer. The program code implemented on a classical computer will be 
a straight connection between the mathematical formulation of quantum me- 
chanics and computational methods. The computational language will include 
formulations such as quantum state, superposition and quantum operator. 

1 Introduction 

The history of quantum computing begins around the 1980s when Richard Feyn- 
man at the First Conference on the Physics of Computation showed that it is not 
effective to simulate an evolution of a quantum system on a classical computer. An 
effective simulation of quantum system has a run-time in polynomial size, i.e. the 
computational time is smaller than a polynomial function of the problem size. The 
conclusion is that relevant simulations of quantum computers will always be larger 
in size than polynomial time. This will bring our attention to the super-polynomial 
time simulations of quantum algorithms; these kinds of simulations have a long run- 
time for large problem. By separating the problems in smaller parts we can avoid the 
long run-time. For examples, there will be super-polynomial time to simulate Shor's 
factoring algorithm on a classical computer. The simulation of quantum algorithms 
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will still be constructive for parts of a larger problem and it will give us a basis for 
comparison between experimental results and theoretical results. The results from 
Shor's algorithm will be verified by multiple the factors from the algorithm outcome 
and hence it is straightforward to check the results from Shor's factoring algorithm 
implemented on a quantum computer. It might be more complicated to check the 
outputs from future algorithms. However, it is possible to show that Shor's algorithm 
will give mathematically correct results, see [6]. But how can we verify that imple- 
mentations of Shor's algorithm on a quantum computer which coincides with their 
mathematical model? A simulation of a quantum algorithm in classical computer 
will give us the possibility to compare the outcome from a quantum computer with 
the output form a physically more stable classical computer. In the development 
of quantum algorithms it will be interesting to check new algorithms on a classical 
computer. This article will show a connection between a future quantum computer 
and today's simulations of quantum computers. The intention is to demonstrate a 
computational simulation of quantum algorithms on a classical computer. A couple 
of examples of various methods for the simulation of quantum algorithms were given 
in [U E]. One of these is the Mathematica package Qdensity which is a simulation 
of a quantum computer developed by B. Julia-Diaz, J.M. Burdis and F. Tabakin. 
Qdensity gives a simulation of a number of well-known quantum algorithms such as 
teleportation, Shor's algorithm, Grover's search and more. An operation in Qden- 
sity is matrix operation and therefore matrixes will represent quantum gates and 
column matrixes will represent a qubit. The package includes the most up-to-date 
essential operations in quantum computing. As examples of these operations we will 
mention the Hadamard, controlled not and Tojfoli gates. In preprint [3] we intro- 
duced a computational language constructed on the basis of quantum mechanics. 
We have decided to implement the well-known Shor's algorithm as an example of a 
quantum algorithm. The aim is to construct a computational language in order to 
present a straightforward connection between Dirac's mathematical formulation of 
quantum mechanics and the program code. Thus, the computational language will 
include quantum mechanical terminology such as quantum operators and quantum 
states. A mathematically described algorithm will have a computational algorithm 
that has a clear mathematical structure in the program code. This simulation will 
be a link from the quantum algorithm described in terms of mathematical concept to 
the implementation of this quantum algorithm on a quantum computer. First we will 
construct a simulation framework in the high-level program language Mathematica. 
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2 The Simulation Framework 



This section will introduce a framework construct for the simulation of quantum 
algorithms in classical computers. We wish to point out that there is a symbolic 
similarity between our framework and the mathematical framework. This framework 
will be a computational dual to Dirac's bra-ket notation. A quantum state in n 
dimensions can be represented by a linear combination of n numbers of basis vectors. 
In the two-dimensional case a quantum state |0) is represented as a superposition of 
two basis vectors, say |0) and |1), (computation basis [ll|5]). In this case a quantum 
state 10) is represented as 



where a and (3 are complex numbers that have the sum of squares equal to one. We 
will introduce some new symbols for the states of this computational basis as follows: 
e[0] = |0) and e[l] = This is the foundation for the structure of the program code. 
For more than one-qubit we will use the computational basis states e[xi, . . . , x„] = 
\xi . . . Xn), where Xj G {0, 1} or by using the more compact notation e[y] = \y), where 
y = x„2° + ■ ■ ■ + xi2^~^. We will, write the state (f) as e[0] = ae[0] + /?e[l], by analogy 
with the equation ([T]). The operator A acts on the state and is often written as 
A\(j)) in the quantum mechanical literature. To match these symbols, we will use 
the computational symbols A|e[0] for this operation. A computational problem is 
to handle different expression with identical meaning. As an example the expression 
e[0,e[l],l] must be equal to e[0, 1, 1] in the code. We can bring in the command 
e[0,e[l],l] := e[0, 1, 1] or the more general e[a__, e[6__], c__] := e[a,b,c] to solve this 
problem. Moreover the program code must be able to handle the linearity of tensor 
product. Let e[.] be vectors and a a scalar. We define the tensor product as 



We add two commands to the program code that will implement this definition of 
the tensor product. The command 



will transform e[a] ® ae[x] ® e[c] to ae[a ® a; ® 6] = ae[a, x, b]. This command is the 
computational dual to the tensor expression in Dirac's notation \a) ^ a\x) <^ \b) = 
a\axb). The other command 




(1) 



Q;(e[f] (S) e[w]) = {ae[v]) (8> e[w] = e[v] (a;e[t(7]) 
{e[vi] + e[v2]) ® e[w] = e[vi] ® e[w] + e[v2] ® e[w] 
e[v] ® ie[wi] + e[w2]) = e[v] e[wi] + e[v] (g) e[w2]. 



(2) 
(3) 
(4) 



e[a__, a_.e[a;__], := ae[a,x,b] 



e[a,^{ae[x] + [3 e[y]), b] := ^ae[a, x, b] + ^/3e[a, y, b] 
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will transform e[a] (S> ^(ae[x] + Pe[y]) (8> e[b] to C,ae[a,x,b] + ^Pe[a,y,b]. Let U be 
an arbitrary one-qubit quantum gate. Then U will transform an arbitrary state e[0] 
which is represented in the computational basis states as e[0] = ae[0] + 6e[l] to the 
state f/| e[(/)] — >• a(cie[0] +C2e[l]) + 6(c3e[0] +C4e[l]), where a,b,Ci G C. Now we 
add a Mathematica gate U to the program code, with the following result f/| e[0] — *• 
ci e[0] + C2 e[l] and e[l] csefO] + C4e[l]. For example, the Hadamard gate H will 
be added in Mathematica as the command if:={e[0] 1/V2{e[0] + e[l]), e[l] — 
1/V2{e[0] - e[l])}. We will define a one-qubit gate Oi as an operator which acts 
on the qubit in position i and leaves the other qubits unchanged. The program 
code must be able to operate with a gate on an arbitrary qubit. Consequently we 
will define an operator in the Mathematica code. The operator will be defined 
Q. _ j(g)«-i (g) f/ (g) j®n-t operator which acts on n-qubits,where / is the one- 

qubit unit operator and U is an arbitrary one-qubit operator. Then operator Oi 
is a function of Oj|e[f] — > e['?/']. Similarly, we will define Ojj as an operator which 
operates as the two-qubits operator on the qubits in positions i, j and leaves the other 
qubits unchanged. Now we have the tools to build quantum circuits. To achieve this, 
we will use a quantum Fourier transform circuit in Shor's factoring algorithm. 

3 An Introduction to Shor's factoring algorithm 

Prime factorizing of an odd number N can be accomplished using Shor's algorithm 
[6]. If is an even integer then we divide it with 2 n-times until 2~'^N becomes an 
odd number. An even N = 2n can easily be found in view of the fact that 2n = 
(mod 2). Let N be the composite of prime factors so that 



where A; > 1 and G Z+. The algorithm will be able to factorize the integer A^. 
We can also assume that is not a prime power, i.e. k > 1 and that there exists 
at least one Pi Pj- A prime power can be found with a known classical method in 
polynomial time. Let us choose a a; G Zn- randomly, where = {1,2,..., A^ — 1, A^}. 
The next step is to use Euclid's algorithm which determines the greatest common 
divisor. If x and A^ are not relatively prime, then we will find a factor by using 
Euclid's algorithm. A factor is equal to gcd(x, A^) if gcd(x, A^) 7^ 1. If x and A^ 
are relatively prime, the task will be to find the order of x in the group Zn. The 
algorithm will search for the smallest r G Z+ so that x'^ = 1 (mod A^); consequently 
r is called order r of x. If 



N=p'i^pr...pT, 



(5) 



x' = 1 (mod A^) 



(6) 
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and r is a even integer, then it is possible to factorize x'^ — 1 as 

- 1 = (xi - + 1). (7) 

The integer N will share at least one factor with (x^ — 1) or (a;5 + 1) since N divides 
— 1. These factors can be calculated with Euclid's algorithm. Let us presume that 
an even r has been found, which gives us the possibility to determine the factors equal 
to gcd(a;5 + 1,A^) and gcd(x5 — 1,A^). There is nevertheless an obvious exclusion 
that X2 = 1 (mod N) due to the definition of order r. However, it may occur that 
x^ = —1 (mod N), i.e. only trivial factors will be found since N is the greatest 
common divisor to X2 + 1 and N. There will be two factors, gcd(x2 + 1,A'") and 
gcd(x2 — 1,N), in case the order r of x is an even number and x^ ^ —1 (mod N). 



4 The Quantum Fourier Transform 

The Quantum Fourier transform is a discrete Fourier transform procedure applied 
to the framework of quantum mechanics, used to realize Shor's algorithm and some 
other really interesting quantum algorithms. The discrete Fourier transform maps 

a vector Vi = ao, ai, . . . , ttn-i to another vector V2 = (Sq, (Si, . . . (3k ■ ■ ■ , (3n-i, where 
a,/? G C. In traditional mathematics the discrete Fourier transform is defined as 
follow: 



n 

j=0 

The quantum Fourier transform (QFT) is analogue to the discrete Fourier transform, 
yet QFT will be represented in computational symbols. The QFT is a function which 
acts on g = log2 n qubits. 

Definition 1. The quantum Fourier transform is a function which maps basis states 
to a linear combination of basis states 

^ n—l 

QFT|eb1 = -^J]e2-^-^/"e[fc]. (9) 

^ jfc=0 

The quantum Fourier transform for an arbitrary state ip is 

n—l ^ n—l 

gFT|5^a,eb1 = ^ ^ a.e^-^^-^/"e[fc]. (10) 
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A decomposition of QFT will be necessary in the task in order to implement 
the circuit on a quantum computer; therefore we will work with the decomposition 
of QFT in our simulation. This decomposition is one of the essential items in the 
creation of a QFT circuit, since it will allow us to construct the circuit with the use 
of quantum gates. Let us introduce the decomposition as 

Lemma 1. 

QFT\ e[j] = ^( e[0] + e[l]) ® ( e[0] + e[l]) ® ■ ■ ■ ® ( e[0] + e[l]). (11) 



We will use the combination of the compact notation e[y] and the notation 
e[a;i, . . . ,Xn] where y = a;„2° + ■ ■ ■ + a;i2""^ and xj G {0, 1}. Hence, we write e[y, 1] 
and e[y, 0] with the meaning e[a;i, . . . , 1] and e[xi, . . . , Xn-i, 0] respectively. The 
proof of this lemma will follow from the fact that e[/,0] G { e[0], e[2], e[4], . . .} and 
e[/,l]G{e[l],e[3],e[5],...}. 

Proof. 

n—l 

QFT\e[j] = — Ve^-^^-'^Z-effc] (12) 

n/2-1 n/2-1 

= ^(V e2-*^'2'/"e[/,0]+ y e2-^^(2'+i)/"e[/,l]) 

^ 1=0 1=0 
n/2-1 n/2-1 

= ^( V e2"*^'2'/"e[/,0] + V e2"'^'2^/V"'^>e[/,l]) 

^ 1=0 1=0 

n/2-1 

= E e'-'^'^/^^^-"^ e[k]) ® (e[0] + e^e[l])) 

v^"" k=0 

n/A-l 

= -^{{Y, e2'^*^''=/("'"')e[fc])®(e[0] + e^e[l])®(e[0]+e^e[l])) 

^ fc=0 

The proof will follow if we continue with this method and decompose the remaining 
qubits in a similarly way. However, the Rotation, the Hadamard and the CNOT 
gates will realize the decomposition in the circuit simulation, illustrated in Figure [H 
Every arbitrary unitary operator may be represented by combinations of single qubit 
gates and CNOT gates, see [5]. QFT can be expressed in single quantum gates as: 

SwapiH,{R,^,_,H,_^{- ■ ■ (i?,,2 ■ ■ ■ R3,2H2iR,,i ■ ■ ■ R2,iH,\ e[u]))))), (13) 
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'q-2, 



Jo / ] H [ ■ ■ iRq-2HRq-1 



H 



Rq-3 " Rq-2 



Swap 



H -Ri 



♦ H 



Figure 1: Circuit for the Quantum Fourier Transform 

where swap is a combination of 3g numbers of CNOT gates. This decomposition 

of QFT requires q operations on the first qubit and q — 1 operations on the second 
qubit, and so on. Hence it follows that the decomposition needs \{q + ^)q H and Rk 
operations. To obtain the right order, we swap the decomposition; thus we need 3g/2 
or 3(g — l)/2 more operations. Altogether, the decomposition of QFT is requires in 
the order of gates e.i. QFT uses 0{q^) elementary operations. 

5 Quantum computation and Shor's factoring al- 
gorithm. 

Shor's algorithm will be divided in four steps, the first, second and last step can 
be executed on a classical computer in polynomial time and we will therefore be 
including these steps as classical algorithm in the simulation. First let us choose an 
arbitrary integer x G Z"^ which will be smaller than the integer N that we want to 
factorize. The second step is to make sure that the chosen integer is not a prime 
factor; otherwise the algorithm has found a factor. It is only in the third step we 
need to implement the algorithms on a quantum computer. In this stage we will use 
the quantum Fourier transform in the algorithm to find the order r of x. Restart the 
algorithm if r is odd or a; 2 = —1 (mod N). In the last stage a classical computer will 
calculate the factors and output gcd(a;2 ± 1, A^). Let us start the algorithm with an 
N and choose n = 2* so that N"^ < n < 2N'^. It will prepare two registers e[0] e[0] 
with q- qubits in the quantum computer. 

1. Let us set up the first register in superposition 

^ n—l 

^j;e[c]e[0]. (14) 

^ c=0 
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2. Let us then compute x'^ (mod A?") in the second register and the computer will 
be in state 



n-1 



—= ^ e[c] q{x'' modA^]. (15) 

^ c=0 

3. Next, we need to compute the QFT on the first register e[c] 

^ n— 1 

QFT\ e[c] = ^ y e='"'^'=/"e[A;]. (16) 



fe=0 

The machine will be in state 

n— 1 n— 1 



^ = - J] J]e2™'=/"e[/t]e[a;'^ modiV]. (17) 



c=0 fc=0 



The algorithm will measure the first register for any A;, but we calculate a measure 
of both registers, since it is the same probability to measure h from the first register 
as there is to measure e[/c, x'^ modA^] from both of the registers. The algorithm will 
find a prime factor only if r is an even integer and the following condition is satisfied 
X2 7^ —1 (mod N\ The probability that a randomly chosen x will satisfy these 
two conditions is at least 9/16. If we choose n > N'^ the quantum algorithm will 
find a useful r with the probability that is at least logiogAf "^^^ algorithm will 
be successful, with high probability if we execute this algorithm for ^^oiogiogAf ^jj^gg 
Now let us measure the first register in the quantum machine for any k in e[k]. The 
order r of x can be found as a denominator of one of the convergents to k/n, where 
the probability to find r depends on the number of qubits. To find the order r, we 
need to apply continued fractions to the k/n 



6 Shor's Algorithm in the Mathematica code 

This chapter will introduce a Mathematica code which implements Shor's algorithm 
on a classical computer. We will follow the Mathematica code evolutions and compare 
this with Shor's algorithm. This comparison will demonstrate a connection between 
the classical computer and the quantum computer. The program code is divided in 
parts with the comments below the code to make it more readable. The program 
will try to find two factors to N, where N is an odd composite of prime factors and 
has at least two different prime factors. 
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N= 3*5; q= rLog[2, N^Jl; 
Do [x = Random [Integer, {2, N-1}]; 
If [GCD[x, N] == 1, SecondStep; QFT; OutPrint 
, Print [ "Chosen x=", x, " a multiplier of " , GCD[x, N] ] , ".";]; 
160Log[Log[N]] ] ^ ^ 
'I 3 IJ 

The prospect for the algorithm to success will increase with numbers of qubits, 
but the execute time will also increase. Of course, this simulation will be a super- 
polynomial time algorithm and large integer will not be able to factorize in reasonable 
time. The algorithm will choose q= \Log2{N'^)'\ so that the algorithm will find a 
factor with large probabihty, i.e. if it selects q to satisfy N'^ < 2^ < 2N'^, the two 
factors will then be found with a probability of at least leolofiogiv' 
The program will choose a random integer x e {2, 3, • • • , A?" — 1}. 

<< NumberTheory^ ContinuedFractions^ 
e[a , a_. e[x ], b ] := ae[a, x, b] 

e[a , (a_. e[x ] e[y ]), b ] := ga&[&, x, b] +Ci8e[a, y, b] 

I v_ := Chop[Expand[v /. {e[x ] ReplacePart [e [x] , e[{x}|[i]]] / . O, i])]] 

J I v_ := Chop [Expand [v / . (e[x ] :■* 

(e[{x}|[i]|, {x} I jj ] /. O /. e[a_, j8_] :-*e [Sequence @@ Take [{x}, i - 1] , 
a. 

Sequence @@ Take [{x}, {i + 1, j-1}], 
/3, 

Sequence @@ Take [{x}, {j + 1, -1}]]))]] 
H := {e[0] « — (e[0] +e[l]), e[l] « — (e[0] -e[l])} 

R := {e[l, e[l, 1]} 

Swap := {e[i_, j_] » e[j, i]} 

We will only use the computational basis states e[a;i, . . . , a;„], where Xj G {0,1}. 
The commands e[a__, Q;_.e[a;__], := ae[a,a;,6] and e[a, ^(a e[a;] + /3e[7/]), 6] := 
^ae[a, a;, h] +^/5e[a, h] will give the program code a connection to the tensor product. 
The command 0_ i_\v- is a one-qubit operator which takes vector v as an attribute 
and operates with O on the qubit in position i. Likewise, the command Ol i_ 
is a two-qubits operator which takes vector v as an input and operates with O on 
the qubit in position i and j. To compute QFT the algorithm requires three gates, 
Hadamard (H), Rotation (R) and Swap (Swap). 
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2 2'J-l 

FirstStep[q_, x_, N_] :=Expand[ ^ 

e [Sequence @@ IntegerDigits[c, 2, q] , Sequence @@ IntegerDigits[0, 2, q] ] ] 

^ 2''-l 



SecondStep 



Secstep[q_, x_, N_] :=Expand[- 



e [Sequence @@IntegerDlglts[c, 2, q] , 
Sequence @@ IntegerDigits [Mod [x°, N] , 2, q] ] ] ; u = Secstep[q, x, N] 

The command First Step prepares the first register in a superposition. Since the 
first step is pointless on a classical computer; consequently it will be excluded in 
the code and we go directly to the second step in the code. Secstep calculates x'^ 
(mod N) in the second register, where q is the number of qubits. 

QFT := (For[i = 1, i ^ q, i+ + , u = Hi | u; For[j = i + 1, j ^ q, + , u = Ri,j | u] ] ; 

For[i = 1, i ^ IntegerPart [ — ] , i + + , u = Swapi^q+i_i | u] ; 

OutQFT = (u / / . a_. e [y ] + b_. e [y ] :-» Together [ (a + b) ] e [y] ) ; 

Probability = List @@ OutQFT / . a_. e[y ] -» {Abs[a]^-, e[y]}; 

Probability = Probability / . {a_, e [y ] } -» a; 

b = {Probabilityll]] } ; 
For [i = 2, i ^ Length [Probability] , i+ + , b = AppendTo[b, b|[i - IJ + ProbabilityliJ ] ] ; 
r = Random [ ] ; 

For[i = 1, i ^ Length [b] , i+ + . If [r <. b[ [i] ] , MeasureQFTStep = i; Break[] ] ] ; 
p= (List @@ OutQFT / . _. e[x ] :-»FroniDigits[ {Sequence @@ Take [{x}, q]}, 2] ) j 



The QFT will act on the state u of the three gates in the following order: 
Hg{Rg,g_iHg_i{- ' ' {Rg,2 ' ' ' ^^3,2 ( i^^, 1 " " " i?2,i^i G [u] ) ) ) ) . The third line in the pro- 
gram code will swap the qubits. All terms with identical computational basis states 
will be collected in the command OutQFT. Probability is a list of the probabilities 
used to measure one of the terms in the register. One of the terms will be randomly 
chosen, taking into consideration of probability to measure the state. The position of 
the chosen term will be saved in MeasureQFTStep. Finally, the hst p of decimal 
numbers is derived from the binary list OutQFT. 
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( r r PlI**®*sureQFTStep]| , , 
OutPrint := ICFD := Denominator [^Convergents [ J J; 

Do[{lf [Mod[x'^™l'^1 , N] == 1 && EvenQ[CFD|[il] &&Mod[x^^^, n] ^N-1, 
Print ["Factors ai = ", GCD[n, + l] , " and", " a.2 = " , 

GCd[n, x^^^-l], " have been found."];]}, {i, 1, Length[CFD] }] ; j 

The randomly chosen value in the register is in p[[MeasureQFTStep]]. In CFD 
the program saves the denominator of convergents p[[MeasureQFTStep]]/2^. From 
this wc can select all even denominators, where x*^^^ = 1 (mod A^) and ^ 
N — 1 (mod N). If any of the denominators satisfies these three conditions, it will 
give us two factors. 



7 conclusion 

In this study we have constructed a computational language for the simulation of 

quantum algorithms and presented the program code for an algorithm. We have also 
demonstrated that every unitary operator has a representation in this computational 
language. An important future challenge is to develop this computational language 
to include the theory of density operator. 
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